//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : solvepsi.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef SOLVEPSI_H
#define SOLVEPSI_H

#include "kernels/solvepsi_kernels.h"

namespace wavescomputer
{

namespace solvepsi
{

template <class T>
inline void phase1(T* bpsi, const T* h0, const T* D0,
            const T* phi, const T* SRCp, const T* SppC,
            const int ld,
            const int Nx1, const int Nx2,
            const int cx, const int cy,
            const T dx1, const T dx2,
            const T pslow)
{
    dimGrid.x = cx / NUM_THREADS_X;
    dimGrid.y = cy / NUM_THREADS_Y;
    dimBlock.x = NUM_THREADS_X;
    dimBlock.y = NUM_THREADS_Y;

    solvepsi::kernel_phase1<T> <<< dimGrid, dimBlock >>>
        (bpsi, h0, D0, 
         phi, SRCp, SppC, 
         ld, Nx1, Nx2,
         dx1, dx2, 
         pslow);
    checkCudaError(cudaGetLastError(), "wavescomputer::solvepsi::kernel_phase1");
}

namespace reflectiveneumann
{

template <class T>
inline void north(T* bpsi, const T* stoN,
           const T* psi_in_N_ext, const T* psi_in_N_cen,
           const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx1 / NUM_THREADS);

    solvepsi::reflectiveneumann::kernel_sp_rn_north<T> 
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoN, 
             psi_in_N_ext, psi_in_N_cen, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::reflectiveneumann::kernel_sp_rn_north");
}

template <class T>
inline void east(T* bpsi, const T* stoE,
          const T* psi_in_E_ext, const T* psi_in_E_cen,
          const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx2 / NUM_THREADS);

    solvepsi::reflectiveneumann::kernel_sp_rn_east<T> 
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoE, 
             psi_in_E_ext, psi_in_E_cen, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::reflectiveneumann::kernel_sp_rn_east");
}

template <class T>
inline void south(T* bpsi, const T* stoS,
           const T* psi_in_S_ext, const T* psi_in_S_cen,
           const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx1 / NUM_THREADS);

    solvepsi::reflectiveneumann::kernel_sp_rn_south<T> 
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoS, 
             psi_in_S_ext, psi_in_S_cen, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::reflectiveneumann::kernel_sp_rn_south");
}

template <class T>
inline void west(T* bpsi, const T* stoW,
          const T* psi_in_W_ext, const T* psi_in_W_cen,
          const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx2 / NUM_THREADS);

    solvepsi::reflectiveneumann::kernel_sp_rn_west<T> 
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoW, 
             psi_in_W_ext, psi_in_W_cen, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::reflectiveneumann::kernel_sp_rn_west");
}

} // namespace reflectiveneumann

namespace sommerfeld
{

template <class T>
inline void north(T* bpsi, const T* stoN,
           const T* psiN, const T* psi_in_N_ext,
           const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx1 / NUM_THREADS);

    solvepsi::sommerfeld::kernel_sp_som_north<T>
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoN, 
             psiN, psi_in_N_ext, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::sommerfeld::kernel_sp_som_north");
}

template <class T>
inline void east(T* bpsi, const T* stoE,
          const T* psiE, const T* psi_in_E_ext,
          const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx2 / NUM_THREADS);

    solvepsi::sommerfeld::kernel_sp_som_east<T>
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoE, 
             psiE, psi_in_E_ext, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::sommerfeld::kernel_sp_som_east");
}

template <class T>
inline void south(T* bpsi, const T* stoS,
           const T* psiS, const T* psi_in_S_ext,
           const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx1 / NUM_THREADS);

    solvepsi::sommerfeld::kernel_sp_som_south<T>
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoS, 
             psiS, psi_in_S_ext, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::sommerfeld::kernel_sp_som_south");
}

template <class T>
inline void west(T* bpsi, const T* stoW,
          const T* psiW, const T* psi_in_W_ext,
          const int ld, const int Nx1, const int Nx2)
{
    int numBlocks = (int)ceil((float)Nx2 / NUM_THREADS);

    solvepsi::sommerfeld::kernel_sp_som_west<T>
         <<< numBlocks, NUM_THREADS >>> 
            (bpsi, stoW, 
             psiW, psi_in_W_ext, 
             ld, Nx1, Nx2);
    checkCudaError(cudaGetLastError(), 
        "wavescomputer::solvepsi::sommerfeld::kernel_sp_som_west");
}

} // namespace sommerfeld

} // namespace solvepsi

} // namespace wavescomputer

#endif // SOLVEPSI_H

